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The freezing transition in a classical three-dimensional system of parallel hard cubes with rounded edges 
is studied by computer simulation and fundamental-measure density functional theory. By switching the 
rounding parameter s from zero to one, one can smoothly interpolate between cubes with sharp edges and 
hard spheres. The equilibrium phase diagram of rounded parallel hard cubes is computed as a function of 
their volume fraction and the rounding parameter s. The second order freezing transition known for oriented 
cubes at s = is found to be persistent up to s = 0.65. The fluid freezes into a simple-cubic crystal which 
exhibits a large vacancy concentration. Upon a further increase of s, the continuous freezing is replaced by 
a first-order transition into either a sheared simple cubic lattice or a deformed face-centered cubic lattice 
with two possible unit cells: body-centered orthorhombic or base-centered monoclinic. In principle, a system 
of parallel cubes could be realized in experiments on colloids using advanced synthesis techniques and a 
combination of external fields. 



I. INTRODUCTION 

Most liquids freeze into a regular period crystalline 
lattice upon a sufficient temperature decrease or pres- 
sure increase. Since this transition is associated with 
a breaking of translational symmetry, it is typically dis- 
continuous (or first-order) .'i' This is in marked contrast to 
two spatial dimensions where no long-range translational 
order existfl and the liquid-solid transition can be con- 
tinuous, for example following the two-stage Kosterlitz- 
Thouless-Halperin-Nelson- Young scenario^. One of the 
few (if not the only) exception to the common finding 
of first-order freezing in three dimensions is a system of 
parallel hard cubes where a disordered liquid freezes con- 
tinuously into a simple cubic (sc) lattice at a volume 
fraction of about 50 percent. This finding was first sug- 
gested for parallel hypercubes in more than three spa- 
tial dimensions by KirkpatriclP by using a second or- 
der virial expansion. Later, Cuesta and coworker^^HH ap- 
plied fundamental- measure density functional theorySEEl 
to parallel hard cubes in three dimensions showing that 
the freezing transition is also continuous in three dimen- 
sions. The continuous nature of the freezing trans ition 
was also confirmed later by computer simulation^^^ 
where the associated criticality was found to be consis- 
tent with the Heisenberg universality class.^i^l 

In this paper, we consider a more general shape of 
hard particles, namely cubes with rounded edges. Our 
motivation to consider rounded cubes is threefold: First, 
it is interesting how persistent the second-order transi- 
tion is with respect to a change of parameters. It is 
known that orientable cubes (i.e. cubes with full ori- 
entational degrees of freedom) freeze via a first-order 
phase transition^^ ^'^ such that additional rotational de- 
grees lead back to the normal picture of freezing. There- 
fore it is interesting to which extent the degree of round- 
ing affects the order of the transition, in particular, since 
the extreme hard-sphere limit d — )• I exhibits the common 
first-order freezing scenario. 



Second, fundamental-m easur e density fu nction al the- 
ory (FMT) was developecPE^ and appliecPE^ further 
in recent years towards hard bodies of arbitrary shape. 
As a version of FMT already exists for the limiting cases, 
parallel hard cubes and hard spheres, rounded cubes con- 
stitute an excellent model system to test the performance 
of fundamental-measure density functional theory. 

Last but not least, it is now possible to fabricate 
micron-sized colloidal particles with almost arbitrary 
shapes^^ ^^. In a recent pioneering work of Rossi et al^^ 
well-defined colloidal cubes were prepared and studied 
in real-space by confocal microscopy. These suspensions 
nicely realize the hard cube model of classical statisti- 
cal mechanics. In the experimental samples, however, 
the cubes typically possess rounded edges, therefore our 
model shape model is closer to these colloids than the 
hard cube. In the experiments, the colloidal cubes are 
not oriented in a parallel fashion. Furthermore, non- 
adsorbing polymers were added to speed up the crys- 
tallization process. Therefore, a first-order freezing tran- 
sition was found in this suspension. Colloidal c ubes c an 
in principle be oriented by external aligning field^^HHl for 
instance by introducing an inner core with two distinct 
non-parallel dipole moments, each of which couples to a 
separate external field. By simultaneously applying two 
non-parallel external fields, which could be external elec- 
tric or magnetic fields or a light field, the orientation of 
the particle described by its two independent axes can 
be fixed. We should note here, that the phase behav- 
ior of a system of parallel monodisperse particles with 
only hard-core interactions at constant packing fraction 
is invariant under scaling of a dimension of the particle 
by a constant factor. Therefore, it is not inconceivable 
that our model of parallel rounded cubes would be real- 
ized in experiments, possibly as stretched rounded cubes. 
Furthermore, our work on parallel, rounded, hard cubes 
provides a good starting point for further studies on col- 
loidal rounded cubes that are not aligned by external 
fields. 

The model for a rounded cube that we use is a sphe- 



rocube, which can be obtained by rounding a cube with 
edge length I. The rounding is done by replacing all edges 
by quartered cylinders of diameter d and the corners by a 
spherical octant such that the curvature is continuous on 
the cube's surface, see Fig. [l] If the diameter d is zero, 
the traditional model of parallel hard cubes is recovered 
while for < d < I we are dealing with truly rounded 
cubes. Finally, in the extreme limit d — > I we obtain the 
hard sphere model where freezing is known to be a first- 
order transition into a face-centered-cubic (fee) lattice.'^ 
By splitting the particle surface into planar, cylindrical 
and spherical parts, we propose a continuous interpola- 
tion between a cube and a sphere which is similar in 
spirit but different in practice to the superball interpola- 
tion used recently by Batten et al.^^^ To abbreviate the 
notation we define a rounding parameter s = d/l, similar 
to Batten et al.'s l/g for the superballs, in the sense that 
both s = 1 and 1/g = 1 denote a sphere, while s = 
and 1/q ^ denotes a cube. The overlaps between two 
superballs can only be detected using an involved nu- 
merical algorithm ,123 which leads to numerical difficulties 
as the superball's shape approaches that of a cubc^^^ In 
contrast, the overlap algorithm for parallel spherocubes 
can be given in a closed and very simple form, as we will 
show in Appendix |Aj Furthermore, the spherocube is a 
very convenient model particle for FMT, since the cur- 
vatures that feature in the theory are constant on the 
spherical, cylindrical and planar sections of the particle's 
surface. Therefore, we have chosen to use the spherocube 
as our model rounded cube instead of the superball. 

We explore the rounded parallel cube model by Monte 
Carlo (MC)^^ and event-driven Molecular Dynamics 
(EDMD)^^ computer simulations and by fundamental 
measure density functional theory of freezin^^Uoi adjusted 
conveniently to the rounded shape. As a simulation 
result, we calculate the equilibrium phase diagram of 
rounded parallel hard cubes as a function of packing frac- 
tion and the degree of rounding embodied in the ratio 
s = d/l. The second order freezing transition known for 
oriented cubes at s = is found to be very persistent 
occurring up to high rounding degrees of about s = 0.65. 
This gives evidence that the second-order freezing transi- 
tion can be seen in experiments on rounded oriented par- 
ticles. The fiuid freezes into a simple-cubic crystal which 
is accompanied by a very large vacancy concentration 
in the emerging solid. At further increasing ratios d//, 
freezing becomes a first-order transition into a sheared sc 
lattice and a deformed fee lattice, where the latter can 
have both orthorhombic (ortho) and monoclinic (clino) 
unit cells. Our simulation data for the continuous transi- 
tion line and for the associated vacancy concentration are 
found to be in qualitative and semi-quantitative agree- 
ment with fundamental-measure density functional the- 
ory. The three novel crystals (sheared sc and the ortho 
and clino variants of deformed fee) can also be confirmed 
experimentally and could be useful for designing new ma- 
terials with novel optical and rheological properties. 

The paper is organized as follows; in section [llj we 
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FIG. 1. (a) A spherocube or rounded cube consists 
of a cube (lightest /gray) surrounded by 6 square prisms 
(darker/yellow), 12 cylinder sections (still darker/light blue) 
and 8 spherical sections (darkest/red). Some sections of the 
outer objects have been removed to show the gray cube, (b) 
Cross section of the spherocube showing the edge length a, 
minimum radius of curvature d/2, and the total width I. 



introduce the rounded cube model in detail. We de- 



scribe the simulation technique in section III while pro- 
viding the background of fundamental-measure density 
functional theory in section |IV[ Results are presented in 
section |V] and we conclude in section IVll 



II. THE MODEL OF ROUNDED PARALLEL 
HARD CUBES 

Our model rounded cube, the spherocube, is a special 
case of the sphero-cuboid introduced by Mulder in the 
context of second order virial theory.'^S' A spherocube can 
be obtained by coating a cube with edge length a with a 
layer of thickness d/2, as shown in Fig. [I] Alternatively, 
it can be constructed by rounding a larger cube with edge 
length I — a+d (dotted rectangle in Fig.[T]D), such that its 
edges obtain a nonzero radius of curvature d/2. We will 
use s = d/l as the shape parameter for the spherocubes. 
The volume of a rounded cube or spherocube is given by 
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In Appendix[Aj we present the overlap algorithm for par- 
allel spherocubes, which is surprisingly simple, especially 
compared to the overlap algorithm for superballs.'^^ 

The thermodynamic state of the system is sometimes 
specified using the pressure P, but mostly using the vol- 
ume fraction or packing fraction rj = v-^cP — v,-cN/V, 
where p is the density, N the number of particles and V 
the volume of the system. The temperature T only serves 
to define the energy unit fc^T, where fc^ is Boltzmann's 
constant. 



HI. COMPUTER SIMULATIONS 

In this section, which consists of four parts, the simu- 
lations that were performed in this work are described. 
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First, we determined candidate crystal structures using 
a recent, but well-tested simulation technique,'^ as sum- 
marized in the first part of this section. After that we 
describe the Monte Carlo (MC) and event-driven Molec- 
ular dynamics (EDMD) techniques. The structural and 
thermodynamic properties we measure are listed in the 
third part and, finally, we describe the methods used to 
determine the phase behavior in Sec. |III D[ 



A. Candidate crystal structures 

We find candidate crystal structures by simulating a 
sin gle un it cell with fully variable box lengths and an- 
g^g j32|33 | ^jjg NPT ensemble, that is the number of 
particles iV, the pressure P and the temperature T are 
held fixed. Using periodic boundary conditions, this unit 
cell is replicated indefinitely to roughly approximate an 
infinite crystal. The final configurations of a number of 
compression series form the unit cells of the potentially 
stable crystal structures. This computationally inexpen- 
sive method has been shown to find all stable crystal 
phases when applied to a system where the phase be- 
havior was already knowrf^ and since then has been em- 
ployed to find candidate structures for a number of novel 
systems t^SEII and also to find close packed structures .^^S 

The variant of the method we use is the following: We 
run a large number of fast compression runs, see Ref. [33 
for details. At the lower pressures, the system samples 
many meta-stable states. As the pressure is quickly in- 
creased, the system gets stuck in one of these states. 
Finally, a nearly perfect crystal is found at very high 
pressure. To distinguish between different crystals we 
use the box shape parameter introduced by De Graaf et 
aZ,'^, which is the average length of the box edges times 
the average area of its faces divided by its volume. The 
states are divided in clusters, such that the box shape 
parameter of each state in a certain cluster deviates less 
than 10^"' from the box shape parameter of at least one 
other state in the cluster. The state which used to repre- 
sent the cluster is the state with the highest density, as 
this will often be the most ordered one. In the remainder 
of this work we will refer to this method for determining 
candidate crystal structures as "unit cell simulations" . 

The small system size allows large fiuctuations in pa- 
rameters such as the density, which ensure that all pos- 
sible states are visited. However, the small system size 
would lead to huge finite size effects, if the results from 
these simulations would be directly used to determine the 
region of stability and other thermodynamic properties 
of the crystals that were found. Rather, this method is 
intended to be used in concert with conventional simu- 
lations (see below), which take crystals formed by repli- 
cating the unit cells obtained from this method as initial 
configurations. 



B. Simulation techniques 

We implemented EDMD simulations for hard rounded 
cubes, which allows us to measure the pressure very effi- 
ciently and to quickly equilibrate the system. Molecular 
dynamics for hard particles is implemented by solving the 
equation of motion exactly. As such, hard particles per- 
form free motion interrupted by instantaneous collisions 
in EDMD. Event driven MD simulations are especially 
fast when collisions can be predicted analytically, such 
as for hard spheres, and also for the rounded cubes stud- 
ied in this work, as described in Appendix \K\ 

Although event driven MD simulations are very fast, 
we found it to be more convenient to use Monte Carlo 
(MC) simulations in the following situations: Many non- 
cubic crystals show a deformation of the unit cell upon 
a change in density or pressure. Moves that change the 
shape of the bojP^ can be easily added to Monte Carlo 
simulation^Sl] account for these deformations. Fur- 
thermore, external potentials, such as the ones required 
for the free energy calculations described further on, can 
easily be accounted for in Monte Carlo simulations, while 
they would make the free motion in between collisions 
too complicated to predict the collisions analytically in 
EDMD simulations. Finally, we want to allow the va- 
cancy concentration to adjust to changes in the density 
or pressure. The simplest way to allow the vacancy con- 
centration to change is to allow the box to change its 
shape in a Monte Carlo simulation. A simulation box 
which starts with Mq = ivi°^ x iV^°^ x ivi°^ unit cells 
with one particle in each unit cell can transform into 
M = Nx X Ny X Nz cells for any integers , Ny and , 
such that M > Mq. The resulting vacancy concentration 
is i^vac = 1 — Mq/M. In practice, we use a N'f^ which is 
50 or 100, such that the vacancy concentration is at low- 
est 1-100/101 ~ 0.01 or 1-50/51 ~ 0.02. Note, that the 
minimal vacancy concentration is an order of magnitude 
larger than the vacancy concentrations in common crys- 
tals (for instance the vaca ncy concentration is of order 
10^* for hard sphere^^HMl) ^ However, the vacancy con- 
centrations in simple cubic crystals of rounded cubes are 
orders of magnitude larger than those of hard spheres, 
as we will see below, which allows us to use this sim- 
ple technique to measure t'vac- However, to keep the run 
time of the simulation limited we have to use consider- 
ably smaller system sizes in the other directions: Nx,Ny 
are either 10 or 15. In the other Monte Carlo simulations, 
we have used a system of approximately 1000 particles 
unless mentioned otherwise. 



C. Thermodynamic and structural properties 

The order parameter m which measures the degree of 
crystallinity was introduced by Groh and Mulder.fiSl It is 
defined using the maximum of the Fourier transformed 
density profile for TOj^, where v = x,y,z denotes a di- 
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rection along one of the Cartesian axis and where the 
density profile is averaged over the two other directions 
before performing the Fourier transform: 



nil. 



maxp^(fc). 

k 



(2) 



The order parameter m is defined by to = (Im^;! + |toj,| + 
|toz|)/3. Because the vacancy concentration and there- 
fore the number of unit cells in each direction can change 
in the variable box length NPT MC simulations (see 
Sec. Ill B I , we do not know the reciprocal lattice vec- 
tor k^x,^, in the i^-direction x^, before hand (However, 
we do know its direction, because the particles do not 
rotate). For this reason, we maximize with respect to k 
in Eq. [2] We can use the resulting length of the wave 
vector ki, to determine the number of unit cells in the v 
direction: N^, = k^L^/27T, where Li, is the length of the 
simulation box in the ly direction. Obtaining Ni, in all 
three directions in this way, we can calculate the vacancy 
concentration i^vac as z/vac = 1 — ^/{NxNyNz). 

We measure the equation of state using one of two 
methods depending on the phase of interest. The pres- 
sure of the fluid phase as a function of density is mea- 
sured in NVT EDMD simulations, while the density as 
a function of the pressure of each of the crystal phases is 
measured using NPT MC simulations. 

We also calculated the mean squared deviation from 
the lattice position (MSD) to compare to the FMT data. 
For system with vacancies, the obvious definition of this 
quantity gives infinity because the particles can easily 
diffuse away from their lattice position by hopping to a 
neighboring, empty lattice site. Instead, we measure the 
MSD from the nearest lattice site. However, we need to 
know the positions of the perfect lattice sites; specifically, 
the shift riyfi of the lattice compared to the lattice which 
has one of its lattice sites in the origin. We use EDMD 
simulations with zero total momentum. This means that, 
when a particle hops a lattice constant aq to, say, the 
left compared to the lattice, the lattice shifts a^/N to 
the right because the center of mass is fixed. As such the 
shift rg drifts with time, and needs to be obtained in the 
simulation before the MSD can be measured. A reliable 
way of determining the shift Tq is to use the phase of 
iTiuiku), which should be equal to exp(ifc^riy,o). 



D. Phase behavior 

We know from earlier work for perfect cubeJ^^ and 
sphereJ^ that the phase diagram contains both first and 
second order melting transitions. For the first order tran- 
sitions we use the highly accurate free energy methods 
that were developed by Frenkel and co-workers.!^ We 
summarize these methods in Sec. IIIID II The second 
or der ph ase transitions were located using finite size scal- 
incjSHIl as described in the section after that. 



1. Free energy calculations 

The Hclmholtz free energy of the fluid is obtained by 
integrating the equation of state from the ideal gas limit: 



fl^.M^^og{pP)+ / dpiP/p-i)/p 



(3) 



where, here and in Appendix [B] the free energy is made 
dimensionless by /* = F/iNkeT) - \og{A^/P), A is the 
(irrelevant) thermal wavelength A = h j \/2Tim^^JtBT , m^c 
is the mass of a particle and h is Planck's constant. 

The free energies of the cry stal p hases are measured 
using the Frenkel-Ladd methocpSl^ in which the free en- 
ergy difference between a crystal and the non-interacting 
Einstein crystal is calculated by thermodynamic integra- 
tion. We have made some modifications to the method 
when applying it to the simple cubic crystal phase to al- 
low for a nonzero vacancy concentration. As these modi- 
fications are similar to the one applied for rotating cubes 
in Ref. [TJ| we leave the details for Appendix [B] As an 
example, the free energy, resulting from the thermody- 
namic integration technique, is shown as a function of 
the vacancy concentration in Fig. [2j Clearly, a finite 
(and quite large) vacancy concentration is found, around 
8%. This is surprising, because the packing fraction for 
this free energy, 77 = 0.53, is rather high compared to the 
critical density rjc — 0.47, as determined using the meth- 
ods described below. We calculated two more free energy 
curves as a function of vacancy concentration and the re- 
sulting vacancy concentrations are shown together which 
the results from the variable box length NVT simula- 
tions and the FMT in Sec.|V] Minimizing the free energy 
with respect to vacancy concentration at every density is 
somewhat cumbersome, so we have used the variable box 
length simulations to determine the vacancy fraction in 
most of this work as described in Sec. IIII CI Once a ref- 
erence free energy f*{po) is known at a certain reference 
density po for each crystal, we integrate over the equation 
of state P/p^ to obtain the free energy at all densities, 
similar to Eq. 

When the free energies of all relevant phases at a cer- 
tain aspect ratio s are known, the coexistence densities 
for a given pair of phases 1 and 2 can be found by solving 
Pi(pi) = P2(p2) and /ii(pi) = P-2{p2), where the pressure 
Pi of phase i is obtained from a fit of the equation of state 
and the chemical potential pi ~ F^/Ni + Pi/ Pi- Solving 
these equations for all possible pairs of phases and finding 
at each density the phase or phase coexistence, which has 
the lowest free energy, the phase diagram can be drawn. 



2. Finite size scaling 

For a second order phase transition, the above method 
can not be used due to very large finite size effects near 
the transition. Near the phase transition, we can use 
finite size scalin^^ to find the properties of the infi- 
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FIG. 2. The free energy of a simple cubic crystal of sphe- 
rocubes with s — 0.6 as a function of its vacancy fraction !/vac 
at 7? = 0.53 for iV ~ 1000 particles. 



nite system. Only the behavior of the order param- 
eter with pressure and system size is required to find 
the location of the transition; no thermodynamic inte- 
gration is required in this case. We use the scaling of 
(m) = N^^^mdPc/P - llN"^) and the Binder cumu- 
lantPc/AT EE 1- (m4)/(m2)2 = U{\PJP - llN""^), where 
Pc is the critical pressure, i>i and 1^2 are finite size scaling 
exponents and rh and U are scaling functions. The expo- 
nents and 1^2 fall into certain universality classes (In 
terms of the critical exponents /3 and v often used in the 
literature, — — /3/3i/ and 1^2 — l/3i^). Groh and Mul- 
dei^^ determined the universality class for the melting 
transition of hard cubes without vacancies to be that of 
the three dimensional classical Heisenberg model, which 
has ~ 0.173 and 1/2 ^ 0.472.1221 We have not found any 
evidence that these exponents are changed when vacan- 
cies are included and, therefore, use these values for the 
exponents also here. 

We used the system sizes N = lO'^, 15^ and 20"^, which 
are large enough that we can neglect corrections to finite 
size scalin^i^ (our large system size is considerably larger 
than that of Ref-fH)). 

A method to determine the critical pressure that does 
not use the values of the scaling exponents consists of 
plotting Upf, which does not depend on system size for 
P = Pc, for a number of system sizes. The point where 
the three curves meet is the critical pressure Pc- We used 
this method for s = and s — 0.5. When N'^^ (to) is plot- 
ted against \Pc/ P~-1\N'^'^ the curves fall on top of a single 
master curve, which confirms the finite size scaling ansatz 
that rh{\Pc/ P — 1\N^^) = {m)N^^ is a universal function 
of the relative deviation from the critical pressure. We 
exploited the universalness of this function to determine 
Pc for the values s = 0.1, 0.2, 0.3, 0.4, 0.6 and 0.65 by 
fitting Pc such that the data for m{\Pcl P -^N"'') fall on 
top of a single master curve for all values of s and system 
sizes considered. The collapse is shown in Fig. [3] Note, 
that the vacancy concentration changes as a function of 
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FIG. 3. The scaled positional order parameter, {m)N"^, of 
systems of A*' hard rounded parallel cubes with s — 0, 0.1 0.2 
0.3 0.4 0.5 0.6 and 0.65 as a function of (1 - P/Pc)N''^ , where 
Pc is the critical pressure (that is a function of s) and ui and 
i/2 are scaling exponents. The system size is A'' = 10'^ for the 
results for s = and s = 0.5, for which = 10^, 15^ and 20^. 
The inset shows {m)N''^ near the critical pressure with a fit 
to the results for s = and N = 10^ subtracted for clarity. 
The pluses and crosses correspond to the solid and dashed 
lines, respectively, of the same color in the main plot. 

the pressure. Consequently, the number of unit cells in a 
certain direction changes discretely in these simulations, 
which is the cause of the noise in Fig. [3j Reassuringly, 
the collapse of the positional order parameters is reason- 
ably good, considering the noise. The inset shows a zoom 
near the critical pressure and a fit to the data for s = 
and N = 10'^ has been subtracted for clarity. Most of 
the data in the inset is indeed scattered around zero for 
a range of pressures near P — Pc, where the exception 
seems to be s = 0.65 (black dashed line). Apparently, 
corrections to scaling are more important for this value 
of s, which is close to the triple point where the simple 
cubic phase is replaced by another crystal phase. There- 
fore, the critical pressure Pc for s = 0.65 is less accurate 
than Pc for the other values for s. 



IV. FUNDAMENTAL-MEASURE DENSITY 
FUNCTIONAL THEORY 

In the framework of density functional theorjl^ the 
equilibrium grand canonical potential is obtained by min- 
imizing the density functional 

m = + y'drp(r) (Kxt(r) - fi) , 

with the intrinsic free energy functional -Flp], the exter- 
nal potential Vcxt(r), the density distribution p(r) and 
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the chemical potential /i. We limit our considerations to 
a one-component system but the theory can easily be gen- 
eralized to multicomponent systems. The functional T[p\ 
naturally separates into two parts, F[p\ — J^id[p] + Jcxc[p] 
with the ideal gas contribution 

-FidM = fesT J drp(r) (log {p{r)A^) - l) , 

where A is the thermal wavelength as defined above. The 
excess free energy J^^xc [p] which contains the information 
of particle interactions is not exactly known, such that 
one has to rely on approximations. 

For hard sphere systems, Rosenfeld's fundam ental 
measure theory (FMT^Pand refined versions 
are currently the most accurate density functional ap- 
proached^. RosenfelcP^ also generalized the FMT 
to arbitrarily convex shaped hard interacting parti- 
cles using the Gauss-Bonnet theorem. His theory 
yields good results only for mildly elongated parti- 
cled^. Recently, Hanse n-Goos and Mecke improved 
these considerationJi^E^ within the so-called extended 
deconvolution fundamental measure theory (edFMT). 
First, we briefly review the results of edFMT and subse- 
quently apply it to spherocubes. 



A. General approach of edFMT 

In the low density limit, we can express the excess free 
energy functional as a second order virial expansion: 



hm J"cxc[p] 

p-i-O 



J I' dr'drp(r')p(r)/(r 



(4) 



Here /(r) — exp ( — (^(r) /ksT)) — 1 is the Mayer function 
with pair interaction potential (f between two particles. 
In the context of hard-body interactions the Mayer func- 
tion simply reads /(r) = —1 for overlapping particles and 
/(r) =0 otherwise. The Mayer function can be deconvo- 
luted into a sum of weight functions Wq, that capture the 
g eomet rical features of a single convex particle, as shown 

/(r) 

-^-^ ^wo<S) W3{r) +wi(S) W2{r) - wi W2(r) (5) 

oo 

-^(-l)JW^P'l(^W^^^'l(r). 

J =2 

Here and in the remainder, we will denote scalar quan- 
tities by X, vector quantities by x and tensors of rank 
j > 2 by Xl^l. The entire set of related scalar, vecto- 
rial and tensorial quantities is referred to as {xq,}. The 
operation (g) of two weight functions is defined by 

i«7 (r) = dr'wa{r') *Wry{r' — r). 



Here, we use the generalized scalar product which is 
meant to be a multiplication for scalar quantities, scalar 



product for vector quantities and trace of the product of 
two matrices for tensors of second order. In general, for 
tensors of jth order we have: 



The geometrical weight functions {wa} are given by 
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In this notation, R(f) is the vector that points along 
the direction f from a certain reference point inside a 
particle to its surface. The principle curvatures ki and 
Kn with corresponding principle directions and v^^ 
are defined at R(f). We also define the mean curva- 
ture H = (ki -|- kii)/2, the Gaussian curvature K = kikh 
and the deviatoric curvature Ak = (kj — kii)/2 for con- 
venience. The normal vector is written as fi. With a^ 
we denote the transpose of the vector or matrix. Con- 
sequently, ab"^ is the dyadic product of vectors a and 
b. 

Finally, we introduce weighted densities as a convolu- 
tion of the density profile with the corresponding weight 
function: 



n„(r) 



(6) 



In terms of weighted densities the low density limit Q 
becomes 

hm J"cxc[/o] = k^T / dr [rtong -I- nin2 - rii • n2 
p^o J 



(7) 



The tensorial weighted density Na should not be con- 
fused with the number of particles N . Here, we truncated 
the tensor expansion after the Jth term and introduced 
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the free parameters Q. In the hmit J ^ oo, the ex- 1. Homogeneous fluid 
act low density hmit is recovered provided Cj = 
The tensorial terms account for the asphericity and van- 
ish only for hard sphere particles. In their original work 
on cdFMT/^ Hansen-Goos and Mecke truncated the in- 
finite sum at J = 2. Subsequently, they showed, for an 
anisotropic fluid of spherocylinders, that renormalization 
by C2 7^ 1 better corrects for the influence of the trunca- 
tion than including higher order tensorial terms without 
renormalizing {i.e. with (2 = 

In edFMT the following ansatz is made for the excess 
free energy 



-^oxcM = A:bT J dr$[{rv(r)}], 

where the free energy density $ solely depends on the set 
of weighted densities Clearly, as p — ^ we should 

recover ([7| . As a result of scaled particle theory and di- 
mensional analysis, the free energy density for truncated 
tensor terms reads 



with 



-no log(l - na) + 



Aiing}) ^ 02({?^2}) 

l-na (1-713)2' 



0l({na}) = "1^2 - ni • n2 



i=2 



02 ({'^2}) = Co (n^N^ 'n2 - 712 n2 • n2 



-7T.2 Tr 



[21n2- 



Tr 



«0 



12U3 



(8) 



(9) 



(10) 



The trace of a matrix Xl^l is denoted as Tr[X[2l]. The 
02 term was introduced by Tarazona^*^ within a dimen- 
sional crossover analysis>^ In the edFMT cq — 3/(167r) 
was chosen. In this way the exact third virial coefficient 
of hard spheres is included. For general convex rotating 
(i.e. non-parallel) particles, the edFMT is exact up to 
the second virial coefficient for the isotropic fluid. 



B. edFMT of parallel hard spherocubes 

We now apply the edFMT to a monocomponent sys- 
tem of parallel hard spherocubes. We chose the model of 
spherocubes as the curvatures of all sections of its surface 
are constant, where we divided the surface of the particle 
into its spherical, cylindrical and flat components, see 
Fig. [1] It is thus sensible to split the weight functions 
into a sum of contributions, each of which is related to 
one of these sections of the surface or, for n^, to the cor- 
responding volume. Since the convolution is a bilinear 
operation, the weighted densities also decompose into a 
sum of terms related to the different components. Due to 
this decomposition we are able to determine the weighted 
densities for each component separately and use a coor- 
dinate system appropriate for its geometry. 



As the most simple case, we first study the mono- 
component homogeneous fluid. It is characterized by a 
constant density profile p(r) = N/V. Consequently, the 
weighted densities {ua} are independent of the position 
vector and read tIq, = pnia, with iria — Jdrwa- The 
integrated scalar weight functions rUa for a — 3, 2, 1, 
represent the volume, surface, mean half widtlPa and 
Euler characteristic of a spherocube. Furthermore, we 
notice that the packing fraction 77 is equal to 713. Thus, 
we can express the weighted densities in dependence of 
the packing fraction as 



77i3 



In the case of the scalar weighted densities we obtain 
Too = 1 

7711 = i (3^ — d) 



777,2 = TTC?^ + 3nd(T 
7)7.3 = Vrc- 
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The n2 vector-type weighted density vanishes for the ho- 
mogeneous fluid, which is a consequence of the Gauss' 
divergence theorem and holds for arbitrarily shaped 
particles As a result, the vector term in the free energy 
density (|8| vanishes. Additionally, the tensor weighted 
density ]v]^' is zero due to the cubic symmetry of the 
particle and the traceless nature of ivj ^ As a result, 
truncation at J = 2 order would leave us only with the 
scalar terms, which does not result in the correct second 
virial coefficient. To improve the theory, we need to ex- 
tend ([5]) at least to the first tensor term that does not 
vanish in the homogeneous fiuid. In the case of sphe- 
rocubes, the first nonzero term contains tensors of fourth 
order, i.e. J = 4. The corresponding generalized product 
reads 



[4] 



8 



9 w 2 



In this way we can determine the free parameter (^4 by 
comparison with the exact second virial coefficient i?2. 
The virial expansion up to the third virial coefficient ^3 
reads 



P 
pk^T 



= 1 



B2 



B 



3 2 



(11) 



with B2 — -ivj-c for spherocubes as can be shown by ele- 
mentary geometrical considerations. On the other hand, 
the compressibility factor is given by the derivation of 
the free energy density with respect to 773: 



pk^T N 97/ (1 - ??)3 



(12) 
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where 



a(,S4j — ^; 

&(C4,Co) = co^^ - a(C4) - 1- 



The Taylor expansion of Eq. ( 12 ) around r/ = is 



l + (a+3)77+ (30+6+6)772 + 0(773), (13) 



where we dropped the (^4 and cq dependences for brevity. 



By comparing the terms of order 77 in ( 12 ) and (111 with 



each other, we obtain 0(^4) = 1. Solving this equation 
for (^4 yields 



C4 



77117712 — 3to3 



[4] 



(14) 



The free parameter (^4 only depends on the shape param- 
eter s. 

In the limit of hard parallel cubes (s = 0), we do 
not recover the equation of state found by Cuesta and 
Marti'nez-Raton,— because for non-spherical shapes the 
third virial coefficient is not exact. We can resolve this 
issue by adjusting the constant Cq, which was introduced 
in ( |Io| ). Equating the third virial coefficient in (11) 
with the one in (13 1 and inserting 0(^4) — 1 
6(^4, Co) = B^v^^ — 9, which is equivalent to 



yields 



Co 



97772/53 



47772 ^''^rc 
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In contrast to i?2/wic ~ 4, the dimensionless third virial 
coefficient B^/v^^ depends on the shape parameter and 
has to be determined for every s. We have numerically 
calculated the third virial coefficient using Monte-Carlo 
integration with an approximate error of ~ 10"''. The 
result i?3(s) smoothly interpolates between the analyt- 
ically known B3 for cubes i?3(0) = 9v'^^ and spheres 
^3(1) = lOv'^^W^ Using the respective third virial coeffi- 
cients of these limiting cases, we recover the FMT equa- 
tions of state of spheres and parallel hard cubes. 



form with lattice constant oq = ((1 — i^vacj^s/v)^ j Pro- 
vided that (1 — ^'vac)'T73 > V^^- With this parametrization 
we can determine the weighted densities according to (|6| . 
As for the homogeneous fluid, we truncate the tensor ex- 
pansion at J = 4. For inhomogeneous density distribu- 
tions, the generalized products of the second and third 
order tensors in general do not vanish. Thus, we need to 
determine the free parameters (2 and (^3. In the special 
case s — 0, the infinite tensor expansion in Eq. ^ with 
J — 00 and (j = (— l)-* can be evaluated analytically for 
the simple cubic crystal, which gives the same result as 
the truncated in Eq. |9] with (C2,C3:C4) = (li~3,4) 
for J = 4. In this work, we require that the truncated 
free energy functional is identical to this analytical free 
energy functional for the simple cubic crystal in the limit 
s — )■ 0. Accordingly, we set (C2,C37C4) = (Ij— 3,C4) for 
spherocubes with finite s, where C4 is given by the value 
for the homogeneous fluid (14 1. Numerical minimization 



of the free energy functional with respect to a and 
yields a continuous freezing transition for s < 0.65. 



V. RESULTS 

A. Crystals and regular close packing 

We found two different types of crystals in the unit cell 
simulations described in Sec. |III A| The crystal found at 
low aspect ratios (for near cubes) resembles a sheared 
version of the simple cubic phase, see Fig. |4]^a). The 
(primitive) lattice vectors at close packing are given by 



ai 




a2 




and a3 



Aa 

' (15) 
The Bra- 

vais lattice of the sheared simple cubic phase (shSC) is 
the rhombohedral lattice. At close packing, the packing 
fraction of the shSC crystal is 



where Aa is given by Aa = d{\ — l/\/2). 



/yshsc = wrc/ {^^ - 3Aa2z + 2Aa3} 



(16) 



2. Simple cubic crystal 

We parametrize the density profile of the crystal by a 
standard Gaussian form given by 

3 

p(Q:,i^vac,r) = (l-!^vac) ^ exp(-a(r - R)^), 

R 

where {R} are the lattice vectors of the prescribed crystal 
structure. We regard i^^ac & [0, 1] as the vacancy concen- 
tration and a € [0, 00) as the Gaussian parameter that 
characterizes the profile. For a simple cubic crystal struc- 
ture the parametrization factorizes and takes a simpler 



The shSC crystal has the highest packing for s < 
0.781133; for higher aspect ratios, a phase similar to 
face centered cubic was encountered, which we called de- 
formed FCC (def FCC), also depicted in Fig.Qa). This 
phase actually consists of two phases between which a 
continuous transition is observed. Decreasing s from 1 
(that is starting with spheres), FCC is deformed such 
that only a base centered monoclinic (BCM) unit cell 
can be recognized (a BCM unit cell can also found in the 
cubic unit cell of FCC). As the aspect ratio is decreased 
beyond 0.825079, a body-centered orthorhombic unit cell 
is found to have the highest packing fraction. Both crys- 
tals can be described using the BCM unit cell (which is 
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FIG. 4. (a) The crystals that were found using the unit cell simulations: a "sheared" version of a simple cubic crystal (shSC) 
and a deformed face-centered cubic crystal (def FCC). For the deformed FCC crystal, the deformed cubic unit cell of FCC 
is indicated in black (yellow and green particles), the body centered orthorhombic (ortho) unit cell in blue (green and blue 
particles) and the base centered monoclinic (clino) unit cell in magenta (red, blue and some of the green particles), (b) The 
packing fractions of the various crystals as a function of the aspect ratio s at close packing. For comparison, the packing 
fraction at close packing for superballs which are mapped onto spherocube with aspect ratio s (see Appendix [c| are included, 
as well as the packing fraction of the simple cubic phase. The inset shows an enlargement of the region where the def FCC 
phase has the highest packing, (c) The lattice vectors of the ortho and clino variants of def FCC. In (b) and the inset of (c) 
the region where the various phases have the highest packing are indicated by the labels "shSC" , "ortho" and "clino" and the 
gray area. 



the most general): 



B. Comparison between FMT and simulations 



ai 



a2 





a2y I , and 

^a2z 




(17) 



where Uii, for i — 2,3 and v — y,z are to be determined 
and the 'base' of the unit cell, on which both particles 
in the unit cell lie, is spanned by ai and a2. The com- 
ponents in general have to be calculated numerically. 
The ones at close packing are plotted in Fig.Qc). 

The packing fraction at close packing for these crystals 
are shown in Fig. |4|^b). For comparison, the close-packed 
simple cubic crystal is also shown. Clearly, the simple 
cubic phase has a lower maximal packing fraction than 
the sheared simple cubic phase for s > and, naively, one 
would think that the sheared simple cubic crystal is more 
stable for all densities. However, we will show below that 
the fluid first transforms into a simple cubic crystal phase 
as the pressure is increased for a large range of s values. 
This shows once again that packing arguments should not 
be used to infer the stable crystal at finite pressures. We 
also included the maximal packing fraction for superballs 
in Fig. Qb) as a function of the s-value of the spherocube 
that has a minimal HausdorfT distance to the superball, 
see Appendix [C) Cube-like superballs show two distinct 
crystal phases at close packing, which are quite similar 
to our shSC and def FCC phases. Clearly, superballs have 
a lower packing fraction at close packing for all values of 
s, which is especially marked at low s, where the flat faces 
of the spherocubes allow a very efficient packing into the 
shSC phase. 



In this section we compare the data from the Monte 
Carlo and event-driven MD simulations to FMT results 
for the simple cubic crystal. In Figs, [sf^a) and (b), we 
show the vacancy concentration as measured in variable 
box length NVT Monte Carlo simulations and compare 
with the results from FMT. Also shown arc two black 
points, which are determined from free energy calcula- 
tions. Reassuringly, these points correspond well to the 
other simulation results. Fig. [sjja) shows the dependence 
on the packing fraction for fixed aspect ratio s = 0. The 
trend of the vacancy concentration from edFMT corre- 
sponds to that of the MC simulations, as does the original 
FMT for hard parallel cubes by Cuesta et aP. However, 
both theories underestimate the vacancy concentration 
considerably. Possibly, fluctuations, which are absent in 
the theory, stabilize crystals with higher vacancy concen- 
trations. The dependence of the vacancy concentration 
on the aspect ratio s is shown in Fig. [sjb). The func- 
tional by Cuesta et ai^ cannot be applied for rounded 
hard cubes with s ^ 0. While the simulation data shows 
only a very weak dependence on s, the theoretical result 
increases dramatically with increasing s, such it actually 
overestimates the data at s > 0.1. Nevertheless, both 
theory and simulations show that the high vacancy con- 
centration is not an artifact of the sharp edges of the 
parallel hard cubes, but remains also when the edges are 
rounded. 

In Fig. [6j the equation of state (P/pfc^r) in the ho- 
mogeneous fluid and in the simple cubic crystal with the 
vacancy concentration of Fig. [5] is shown. Again, simula- 
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FIG. 5. (a) The vacancy concentration i^vac as a function of 
packing fraction rj as obtained from variable box edge length 
NVT simulations for parallel hard cubes, that is with s = 0, 
from Cuesta et al 's FW^ for parallel hard cubes and the 
edFMT of this work for rounded cubes with s = 10~*. (b) 
The same as (a), but now s varies and r] is fixed to 0.5 or 
0.55. The black points is determined by minimizing the free 
energy with respect to the vacancy concentration for s = 
and T) — 0.47 in (a) and for s — 0.6 and rj = 0.5 in (b), see 
Fig-i 



tion results are compared to FMT results obtained using 
the functional from this work and from Ref . i5i The agree- 
ment with the simulation data is reasonable for s ~ 0, 
while the previous functionaF" for parallel hard cubes de- 
scribes the simulation data the best. As s increases, the 
agreement between theory and simulations deteriorates 
somewhat, and for large values of s improves again, such 
that the difference in the compressibility factor resulting 
from the theory and the simulation is of the order of 1 
for s > 0.1. 

The root mean squared deviation from the nearest lat- 
tice site (RMSD) as measured in event-driven molecular 
dynamics simulations (EDMD) is compared to FMT re- 
sults in Figs. ^B,) and (b). Again, the first of these fig- 
ures shows the t] dependence at s = 0, while the second 
shows the s dependence at fixed rj. The RMSD is often 
compared with the Lindemann criteriorP^ for first order 
phase transitions, which says that the crystal starts to 
melts when the RMSD is around 10%-20% of the lat- 



FIG. 6. (a) The compressibility factor P/pknT as a function 
of packing fraction rj as obtained from event-driven MD sim- 
ulations for parallel hard cubes (s = 0), from Cuesta et aVs 
for parallel hard cubes and the edFMT of this work for 
rounded cubes with s = lO"*. The FMT and edFMT results 
for the fluid are exactly equal, (b) The same as (a), but now 
s varies and is fixed to 0.3, where the fiuid is found, or 0.5, 
where the simple cubic crystal is stable. 



tice constant. As expected for second order phase transi- 
tions, the RMSD of the simple cubic crystal of (rounded) 
cubes is higher than the Lindemann parameter at the 
transition; it is in fact two to four times as high. The 
theoretical RMSD results at s = [Fig. ^a)] again show 
the correct trend, but both theories under- estimate the 
simulation results, as was the case with the vacancy con- 
centration. The result from cell theory (for zero vacan- 
cies) is also indicated by the think line. Interestingly, 
the simulation results have a substantially different slope 
than cell theory even when approaching close packing, 
the MSD from simulations is approximately 1.6 times the 
cell theory result. For comparison, the mean squared dis- 
placement measured in simulation of hard spheres is ap- 
proximately 1.098(4|^ times the cell theory result. The 
dependence of the RMSD on the s of FMT in Fig. ^h) 
is qualitatively very similar to the simulation results: At 
T] — 0.55 the RMSD decreases monotonically, while the 
RMSD for 77 — 0.55 shows a strong decrease with increas- 
ing s for small s followed by a small increase when s is 
increased beyond a certain value s ~ 0.55. 
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FIG. 7. The dimensionless root mean squared deviation 
(RMSD) from the nearest lattice site divided by the lattice 
constant, Ar/ao, as obtained from variable box edge length 
NVT simulations. Results from our edFMT and the older 
FMT are also shown, (a) The RMSD as a function of packing 
fraction rj for parallel hard cubes, that is, for s — 0. The thick 
black line denotes the cell theory result without vacancies, (b) 
The RMSD for varying s at r? = 0.5 and rj = 0.55. 

C. Phase diagram 

The phase diagram of parallel hard rounded cubes 
is shown in Fig. [sj The fluid to simple cubic (SC) 
crystal second order transitions for hard parallel cubes 
s = from this work have a critical packing fraction, 
rjc = 0.469(3) from simulations and rjc — 0.3325 from 
FMT, which should be compared to the earlier simula- 
tiorfi^ critical density, r]^ — 0.53(1), and the result from 
the original FMT,E1 77 = 0.3143. Our simulation results 
differ from the previous work because it was assumed 
that the SC crystal had zero vacancies in Ref. [121 In con- 
trast, we find an extremely high vacancy concentration of 
13% at coexistence. The FMT vacancy concentration of 
the earlier FMT^'was 30%, while our edFMT gives 23%. 
Note, that the critical densities differ for the two theo- 
ries and the simulations, which explains the reversal of 
the trends compared to the results at fixed packing frac- 
tion. Our FMT describes our simulation results slightly 
better than the earlier FMT,^ as far as the critical density 
and the vacancy concentration at rjc are concerned. Con- 
versely, the inclusion of vacancies has brought the critical 



density from simulations closer to the FMT results. 

For hard parallel cubes with a finite rounding {i.e. 
s 7^ 0), the transition from the fluid to the SC crystal 
phase is also second order both for the simulations and 
for FMT. For the simulations, this is indicated by the 
critical scaling whose exponents belong to the Heisen- 
berg universality class.^ The reasonable agreement be- 
tween the simulation results and the FMT at s = is 
even slightly improved for s 7^ 0, see Fig. [8] The phases 
that were found using the simulations of single unit cells 
all have their separate area of stability in the phase dia- 
gram. At low s, the SC is stable at low densities, while 
the sheared variant (shSC) is stable at high densities. 
The transition from SC to shSC seems to become more 
weakly first order as s approaches zero, and simultane- 
ously the vacancy concentration decreases. For s < 0.5, 
we did not use free energy calculations, because the free 
energies of SC and shSC were very close making it hard to 
find the transition. Instead, we used direct simulations, 
which always lead to a pressure at which the difference in 
chemical potential was very small. As s is increased, the 
SC-shSC transition goes down in density and at some 
point the shSC coexist directly with the fluid. Finally, 
the deformed FCC phase is stable for sphere-like parti- 
cles (s > 0.8). As mentioned in Sec. |V A[ the def FCC 



phase actually consists of two crystals, one with a base- 
centered monoclinic (clino) unit cell and an other with a 
body-centered orthorhombic (ortho) unit cell. Our mo- 
tivation for investigating the clino to ortho transition in 
this system, is that pyroxene, the second most abundant 
mineral in the earth's mantle, also has clino and ortho 
forms.'^ As a result, the clino-ortho transition of pyrox- 
ene is a topic of great interest in geology. In our case, 
the ortho-unit cell needs only be slightly deformed to 
form the clino unit cell: the angle between two of the 
orthorhombic lattice vectors is changed to slightly to a 
little more or less than 90 degrees; the difference is at 
maximum 2.22 degrees for s — 0.884 at close packing. 
The packing fractions at the transition from clino to or- 
tho and the transition from shSC to ortho were obtained 
by direct simulations, as shown in Fig. |8] The shSC- 
ortho phase transition is more strongly first order than 
the ortho-clino transition, which enabled us to calculate 
the free energy of the shSC and ortho phases separately 
and, reassuringly, the free energy difference between the 
two phases at the shSC-ortho transition, which was ob- 
tained in direct simulations, is smaller than the statistical 



VI. CONCLUSIONS 

We studied a system of parallel rounded cubes (sphe- 
rocubes) using fundamental measure theory and simula- 
tion with a special emphasis to the second order freezing 
into the simple cubic phase. We developedthe fundamen- 
tal measure theory starting from edFMT^l^IlSl expanding 
up to fourth order tensor terms and renormalizing the 
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FIG. 8. (a) The phase diagram of parallel rounded cubes in the rj-s representation, where ri = VrcN/V is the packing fraction 
with Vrc the volume of a rounded cube and s — d/l is the rounding parameter (see Fig. [T]). A cube has s = and a sphere 
s = 1. Shown are the areas of stability of the deformed fee phase of near spheres (def FCC), the sheared cubic crystal (shSC), 
the simple cubic crystal (SC) and the fluid phase in white. The forbidden region above the close packing density is shown in 
dark gray and coexistence areas in lighter gray (coexistence lines are vertical). The filled symbols (MC simulations) and the 
thick line (FMT) denote second order phase transitions, while the empty symbols denote first order phase transitions from 
simulations, (b) An enlargement of the large s region of the phase diagram: the def FCC phase is actually seen to have a 
body-centered orthorhombic variant (ortho) and a base-centered monoclinic variant (clino), as depicted in (c). 



third and fourth order terms. When we apply the theory 
to the simple cubic phase, we find that the freezing is 
second order, not just for perfect cubes with shape pa- 
rameter s — 0, but also when we introduce a degree of 
rounding s up to s = 0.65. This finding is confirmed by fi- 
nite size scaling techniques using Monte Carlo (MC) sim- 
ulations. Furthermore, we find both in theory and simu- 
lations an unusually high vacancy concentration, namely 
13%, that is, about twice as high as for rotating per- 
fect cubeJ^ and fo ur ord ers of magnitude higher than 
that of hard spheres .1^^!^ The very high vacancy concen- 
tration and the simple overlap criterion make this system 
an ideal system to study vacancies (the number of vacan- 
cies can always be decreased by increasing the density if 
so required). 

When comparing the theory to the simulations, we find 
good qualitative agreement; exceptions are the depen- 
dence of the vacancy concentration and the pressure on 
the shape parameter s which show the incorrect trend 
for FMT. Quantitative differences between the FMT and 
simulation results are found for most quantities. This 
is most likely caused by anomalous higher virial coef- 
ficients for this system (parallel cubes for instance have 
negative sixth and seventh virial coefficient^^ which are 
not reproduced by the theory. Nevertheless, the most 
important property, the packing fraction at freezing is 
predicted quite well by FMT over the whole range of s 
where the simple cubic crystal is stable. 

Finally, we completed the phase diagram of rounded 
cubes by investigating the possibility of other crystal 
phases in direct simulations and by performing free en- 
ergy calculations using the results obtained in these sim- 



ulations. The phase diagram of parallel rounded cubes is 
surprisingly rich considering the simplicity of the model: 
it contains, apart from the simple cubic crystal phase, 
three more crystal phases. For low values of s and high 
densities a sheared variant of the simple cubic (shSC) 
is found, while at low densities the simple cubic crystal 
is stable. For higher values of s, first the simple cubic 
crystal and later also the shSC phase disappears to be 
replaced by a body-centered orthorhombic crystal which 
is essentially a slightly deformed face-centered cubic crys- 
tal. Finally, a base-centered monoclinic crystal is found 
for values of s near one, that is, for near spheres. The 
resulting symmetry change is interesting because a sim- 
ilar transition is found for an abundant mineral in the 
earth's mantle.^ We expect that most, if not all, of the 
crystal phases we observed can also be found for rotating 
rounded cubes, which could be verified in experiments on 
colloidal rounded cubes. 
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Appendix A: Overlap algorithm and collision 
prediction 

The criterion for overlap between two (co-aligned) 
spherocubes is surprisingly simple: Two spherocubes 
overlap when shortest distance, Ar^j, between any point 
on the surface of particle i and any point on j is smaller 
than d. The shortest distance can be calculated in the 
following two steps: 



spring, such that the external coupling potential reads: 



(Al) 



J sign(rj> - ri^u)hj,u h^^^ > 
Ar„-,^ = < ^ ^, . (A2) 

otherwise 



the norm of the vector with components Ar^^i, for v — 
x,y,z is the shortest distance Ar^. The simplicity of 
the overlap criterion for spherocubes is a large advan- 
tage compared to superballs for which it can not be as- 
certained using analytical means whether two particles 
overlap or not.^^ Note, that the overlap criterion for sphe- 
rocubes becomes more complicated when the particles 
are not aligned with the Cartesian axes. 

Collisions can be analytically predicted for parallel 
hard rounded cubes as follows: The surface of a rounded 
cube consists of sections of axis-aligned cylinders, planes 
and spheres, for which collisions can be easily calculated 
by the collision detection algorithms for the correspond- 
ing one, two and three dimensional hyper-spheres.^^ 
Specifically, a collision test for two n-dimensional hyper- 
spheres is required when exactly n components of hij 



[the vector with components &ij>, see Eq.( Al) ] are 



nonzero (at least one of components of hij is non-zero 
initially in an overlap-free configuration). The time at 
which a certain component of hij becomes nonzero can 



easily calculated from Eq. (All. All such times are de- 



termined, sorted and inserted in a list, whose subsequent 
elements define the time intervals at which certain parts 
of the surface of one of the particles might collide with a 
part of the other particle's surface. For each such a time 
interval, the corresponding hypersphere collision check 
is performed and the shortest of the resulting times is 
the time of collision of the two spherocubes. The rest 
of the algorithm is the same as the optimized algorithm 
for hard sphereSj'^SJ in which the collisions are stored, to- 
gether with others event (such as measurements) in a 
binary tree leading to a theoretical iVlog(A)scaling of 
the computational effort for a fixed run time.l^ 



Appendix B: The Prenkel-Ladd method for 
crystals with vacancies 

In the original Frenkel-LadcF^ approach, each particle 
is coupled to its ideal lattice position with a harmonic 



N 



/3C/har(r'^;A) = A^(r,-ro,OV^' 



(Bl) 



where r.^ denotes the position of particle i and ro.i the 
lattice site of particle i and /3 — l/ksT. If the value of A 
is high enough or if the lattice positions are far enough 
apart, the particles do not interact and the free energy of 
the system is given by the known analytical free energy 
of the non-interacting Einstein crystal.EH Therefore, the 
coupling constant A can be used to switch between an 
ideal Einstein crystal for high A and the unperturbed 
crystal for A = 0. The free energy of the crystal for 
A = can then be found by integrating over A: 



r(A,F,T) = /^i„,t(iv,y,r) 



dA 



(B2) 



where (9/V9A) = (t/har(r^; A))/(AiV). For A = 0, (U) 
diverges as the center of mass of the system diffuses as a 
whole, taking the particles ever further away from their 
lattice position. To overcome this problem, the center 
of mass is fixed which results in additional (small) terms 
in the free energy of the n on-interacting system, which 
can be found in Refs.'^^'^. The value of Amax required 
to obtain a non-interacting Einstein crystal depends on 
the lattice spacing, such that free energy calculations 
over a wide range of densities require constant tuning 
of Amax- However, the same value for A^ax can be used 
for every density if the inter-particle potential is replaced 
by a purely repulsive finite potential whose interaction 
strength is slowly decreased from essentially infinite to 
zero (where an essentially infinite interaction strength 
implies that no overlap is found during the simulation). 
The free energy difference between the interacting crys- 
tal for A = Amax and the non-interacting Einstein crystal 
is then obtained by integrating over the strength of in- 
teraction 7 of the inter-particle potential, see Refs.'S^M] 
for details. The soft interaction between two particles. 



in this case, reads 7(1 — O.QAr^j/l), see Eq. (A2|. We 



have used this method for all crystal phase witn excep- 
tion of the simple cubic crystal phase which had a large 
concentration of vacancies. 

When the crystal has a nonzero vacancy concentration, 
the Frenkel-Ladd methocPi' requires some modifications. 
First of all, the particles are no longer associated with 
a single lattice site as they can hop to a different site 
when it is empty. Therefore, the harmonic potential can 
no longer be used. Instead, we use the periodic external 
potential proposed by Groh and Mulder,'i2l which reads 



/3[/per(r^) = A 



N 

E E 

i—l iy—x,y,z 



1 - cos{kj,ri^j,), 



(B3) 



where fc^ ~ 2ttN^/L^. The number of unit cells N^, in 
direction v and the length of the edge of the box 
are adjusted to tune the density and vacancy concen- 
tration. Furthermore, we promote hopping of a particle 
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from a filled lattice site to an empty one also at large A 
by performing moves of exactly one (cubic) lattice vec- 
tor, which allows the distribution of vacancies over the 
lattice to equilibrate. Note, that the center of mass of the 
system is already fixed by the external potential, so no 
additional terms due to the fixing of the center of mass 
arise in the free energy. However, additional moves which 
translate the whole system homogeneously are required 
to efficiently equilibrate the center of mass. Finally, the 
ideal Einstein crystal free energy itself is modified be- 
cause of the modified inter-particle potential. Addition- 
ally, the combinatorial free energy of choosing N filled 
lattice sites out of a total of M lattice sites needs to be 
included in the free energy. The total free energy of the 
non-interacting Einstein crystal reads 



/Einst — 



Ml 



{M - Ny.m 



31nzi(A„,ax) (B4) 



where 2i(A) = /o(A) exp(— A) ao/A is the factor in the 
partition sum which results from the integration of the 
degrees of freedom of a single particle in one direction, 
and Iq is the zeroth modified Bessel function of the first 
kind. If A goes to infinity, the external potential (B3) 
can approximated by a harmonic potential, and simulta- 
neously 2;i(A) approaches ^/\J2tt a^/ k as Amax oo. It 
is this expression for zi(Aniax)j which we use in practice, 
because the approximation has a negligible effect for the 



direction 


dsh 


drc 


(1,0,0) 

(1.1.0) /V2 

(1.1.1) /^/3 


r 


1/2 

a/V2 + d/2 
VSa/2 + d/2 



value of Amax we used (An 



4000). 



TABLE I. The distance between the center of a particle and 
its surface in the indicated directions for superballs (dab) and 
rounded cubes or spherocubes (d^c)- 



the normals to the respective surfaces at these points are 
equal. 

We only need to consider the distance from the center 
to the surface in this case in three different symmetry 
directions. These distances are listed in Tbl. HI To cal- 
culate the Hausdorff distance, we only need to maximize 
over the directions n given in Tbl. |T} 

d(sh[r,q],Tc[l,s]) = rRa.x\dsh{r,q) - d,c{l,s)\, (C4) 

n 

where we use sb[r, q] and tc[1, s] to denote a superball and 
rounded cube, respectively, with a certain aspect ratio (s 
or q) and linear size {I or r). Minimizing d(sh[r, g], rc[/, s]) 
for a certain value of s with respect to r/l and q, we ob- 
tain the superball that best fits a spherocube with aspect 
ratio s. 
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